%% Description: The SVJEJ2 model
%%
clear;
load('optiondata.mat');
Price=100;Times=4/365;Rate=0.01;DividendYield=0;
%% fit the SVJEJ2 model
options=optimset('Display','iter','MaxFunEvals',10000);
v0=[0.5 0.2 2 0.8 -0.8 2 -0.05 0.03 0.5 -1 0.5 0.5];
lb=[0 0 0 0 -1 0 -2 0 0 -10 0 0];
ub=[1 1 7 5 1 15 2 1 1 10 2 2];
fn1=@(v)rmse_SVJEJ2(optiondata,Rate,Times,Price,DividendYield,v);
[optpar1,minn1]=fmincon(fn1,v0,[],[],[],[],lb,ub,[],options);
%% calculate implied volatilities for the SVJEJ2 model
V0=optpar1(1);ThetaV=optpar1(2);Kappa=optpar1(3);SigmaV=optpar1(4);RhoSV=optpar1(5);JumpFreq=optpar1(6);
mj=optpar1(7);JumpVol=optpar1(8);prob=optpar1(9);mlow=optpar1(10);slow=optpar1(11);shigh=optpar1(12);

Strike=optiondata{:,1};N=length(Strike);
CallPrice_SVJEJ2=optbySVJEJ2(Strike,Times,Price,Rate,DividendYield,V0,ThetaV,Kappa,SigmaV,RhoSV,JumpFreq,mj,JumpVol,prob,mlow,slow,shigh);

IV_SVJEJ2=zeros(N,1);
for j=1:N
    IV_SVJEJ2(j)=blsimpv(Price,Strike(j),Rate,Times,CallPrice_SVJEJ2(j),[],DividendYield);
end
%% fit the DJKS (2019) model
v0=[0.5 0.2 2 0.8 -0.8 2 -0.05 0.03 0.02];
lb=[0 0 0 0 -1 0 -2 0 0];
ub=[1 1 7 5 1 15 2 1 2];
fn2=@(v)rmse_SVJEJ(optiondata,Rate,Times,Price,DividendYield,v);
[optpar2,minn2]=fmincon(fn2,v0,[],[],[],[],lb,ub,[],options);
%% calculate the implied volatilities for the DJKS (2019) model
V0_=optpar2(1);ThetaV_=optpar2(2);Kappa_=optpar2(3);SigmaV_=optpar2(4);RhoSV_=optpar2(5);JumpFreq_=optpar2(6);
mj_=optpar2(7);JumpVol_=optpar2(8);shigh_=optpar2(9);

CallPrice_SVJEJ = optbySVJEJ2(Strike,Times,Price,Rate,DividendYield,V0_,ThetaV_,Kappa_,SigmaV_,RhoSV_,JumpFreq_,mj_,JumpVol_,0,0,0,shigh_);
IV_SVJEJ=zeros(N,1);
for j=1:N
    IV_SVJEJ(j)=blsimpv(Price,Strike(j),Rate,Times,CallPrice_SVJEJ(j),[],DividendYield);
end
%% Plot implied volatilities similar to Figure 8 
plot(Strike,optiondata{:,2},'o',Strike,IV_SVJEJ2,'k',Strike,IV_SVJEJ,'k--')
legend('observed implied volatility','SVJEJ2','SVJEJ')
%% Bimodal PDF generated by the SVJEJ2 similar to Figure 9
xpdf=-0.2:0.001:0.2;
ypdf=DensSVJEJ2(xpdf,Price,Rate,DividendYield,Times,V0,ThetaV,Kappa,SigmaV,RhoSV,JumpFreq,mj,JumpVol,prob,mlow,slow,shigh);
plot(xpdf,ypdf,'k')
